function [bv,sebv,R2v,R2vadj,fit,e] = robustOLS(y,X,k,se_scheme)

%scheme: 1: Newey West (1987), 2: Hansen-Hodrick (1980)

T=size(y,1);
K=size(X,2);

bv=X\y;
fit = X*bv;
e=y-X*bv;

sebv=NaN(K,1);

R2v=1-var(e)/var(y);
R2vadj=1-var(e)/var(y)*(T-1)/(T-K);

d=1/T*(X'*X);

%1: NW and 2: HH
if se_scheme == 1 || se_scheme ==2;

w=e*ones(1,K).*X;
C0=1/T*(w'*w);
S=C0;

for j=1:k;

    Cj=1/T*w(j+1:end,:)'*w(1:end-j,:);
    
    if se_scheme ==1;
    S=S+(k+1-j)/(k+1)*(Cj+Cj');
    elseif se_scheme ==2;
    S=S+(Cj+Cj');   
    end;
    
end;

inv_d = d\eye(K);

varb=1/T*inv_d*S*inv_d;

sebv=sqrt(diag(varb));

end;

